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FOREWORD 


This report is a technical summary of the progress made by the 
Electrical Engineering Department, Auburn University, toward fulfill- 
ment of Contract NAS8-20104 granted to Auburn Research Foundation, 
Auburn, Alabama. This contract was awarded April 6, 1965, by the 
George C. Marshall Space Flight Center, National Aeronautics and 
Space Administration, Huntsville, Alabama • 



SUMMARY 


A numerical integration scheme for solving the four parameter 
vector differential equation is derived and investigated in this report. 
The results obtained can be applied to a large class of numerical 
integration schemes, since this class can be shown to be equivalent 
to the derived scheme. 

Bounds for the truncation errors and roundoff errors generated 
by the digital computer in computing the four parameters using the 
derived scheme are developed. Study shows that the resulting error 
bounds are useful in the determination of an optimal integration 
scheme and sensor sample rate for a particular mission using a 
given computer. 
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I . INTRODUCTION 


Navigation is that branch of art or science of directing the 
course of vehicles. It involves the knowledge of present position, and 
the direction and magnitude of motion with respect to other reference 
points* Host navigation systems depend upon some external aid in ob- 
taining this information, while inertial navigation systems are capa- 
ble of deducing all this information from on-board measurements in 
self-contained system. These on-board measurements are obtained 
means of sensors, such as accelerometers and angular rate gyroscopes 
mounted to the vehicle. There are two methods in mounting these 
sensing devices: the stabilized platform method and the strapdown 

method. 

In the stabilized platform method, the sensors of the inertial 
navigation system are mounted on a stable platform. The platform is 
kept inertially aligned with a predetermined set of inertial axes by 
suspending in a system of gimbals. Therefore the resulting measure- 
ments are in the inertial coordinate system. 

In the strapdown method, the sensors of the inertial navigation 
systems are rigidly fixed to the vehicle and hence the resulting mea- 
surements are in the vehicle coordinate system. Since navigation 
equations are usually solved in the inertial coordinate system, it is 
necessary to generate a coordinate transformation matrix that can in 
turn be used to transform the measured acceleration vector in the 
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vehicle coordinate system to the inertial coordinate system. The 
coordinate transformation matrix is generated by an on-board digital 
computer. This computer utilizes angular rates obtained from the sen- 
sors, which are mounted to the vehicle, to compute the transformation 
matrix. 

The coordinate transformation matrix, C, relating the inertial 
coordinate system to the vehicle coordinate system is given by [6] 

v x - CV^ (I-D 

where V_ v is a column vector with components measured in the vehicle 
coordinate system. 

is the same vector with components measured in the inertial 
coordinate system, and 

C is the square matrix of direction cosines of the inertial 
axes relative to the vehicle axes. 

There are three basic methods of representing the transformation 
matrix C. These methods are 

1. Direction cosine. 

2. Euler angles (three and four angle methods). 

3. Four parameter methods (Euler parameters, quaternions, and 
the Cayley-Klein parameters) . 

In each case the transformation matrix can be computed using a 
set of first order differential equations which require as inputs the 



3 


measured angular rates about the three vehicle coordinate axes. The 
four parameter method will be considered in this study since it has 
fewer computer operations required for its implementation than the di- 
rection cosine method and has no singular point as does the three 
Euler angles method. 

As shown in 11] , there are three different methods (Euler T s theo- 
rem, quaternions, and Cayley-Klein) in deriving. the same coordinate 
transformation matrix as expressed by the four parameters. These meth- 
ods also lead to the same set of first order differential equations re- 
lating the vehicle body angular rates to the time rate of change of the 
four parameters. 

The four parameters may be defined by the application of Euler 1 s 
theorem, which states that any real rotation may be expressed as a ro- 
tation through some angle, about some fixed axis, as 

22 Cos y/2 
e? - Cos y Sin y/2 

(1-2) 

= Cos 3 Sin y/2 
64 - Cos y Sin y/2 

where y is the angle of rotation and a, 0 and Y are direction 

angles between the rotation axis and x, y and z axes of the inertial 
coordinate system. The transformation matrix relating the initial 
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coordinate system and the vehicle coordinate system in terms of the 
four parameters is 

“l"*2"*3 + *4 2< WV2> 2(e l e 3 + V4 ) 

^ ^ e l e 2”** e 3 e 4^ e i~ e 2"*" e 3” e 4 ^^ e 2 e 3~ e l e 4^ (I - 3) 

2 ( e 2 e 4" e i e 2) 2 ( e 2 e 3 +e l e ^ e l +e 2~ e 3~ e 4 

The time rate of change of the four parameters in terms of the ■body- 
rates is 



®1 = 2 ( “^ e 2 " ^ y e 3 " * x e 4 } 
k 2 = | (+ ^ e l " ^x e 3 + 

e 3 = i(+‘^ y e 1 + l x e 2 - 

t *] . * • • 

e. - — (+4> e 1 - <{> e„ + <}> e -) 
4 2 x 1 7 2 z 3 


(1-4) 


♦ % « 

where <j> x , <j> and <j>^ are the measured angular rates of the vehicle with 
respect to the inertial coordinate system . 

Now equation (1-4) is to be solved by various numerical integra- 
tion techniques using an on-board computer to update the four parame- 
ters, which in turn are used to compute the coordinate transformation 
matrix. In order to select an optimal integration scheme, to determine 
the computer sizing and to evaluate the performance of the system re- 
quirements, it is necessary to determine the error introduced by the 
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computational process. The object of this study is to investigate the 
computational errors introduced in computing the numerical solution of 
equation (1-4) using a digital computer. 

The main body of this study is divided into five chapters and five 
appendices. The layout of subsequent material is as follows: 

Chapter IX derives the exact solutions for the four parameters 
\tfhen the angular rates are proportional to each other, A four parame- 
ter algorithm is then presented for error analysis purposes. 

Chapter III analyzes both the truncation and the roundoff errors 
introduced in the digital computation of the four parameters using the 
algorithm developed in Chapter II, Roundoff error bounds for the basic 
arithmetic operations are discussed. Techniques for determining the 
propagated truncation errors and accumulated roundoff errors are de- 
scribed. 

Chapter IV presents the results of two selected examples. 

Finally, Chapter V embodies the conclusions and recommendations. 
Appendix A describes the application of the Peano-Baker method of 
successive approximation. Appendix B discusses the vector and matrix 
norms. Appendix C proves that $(m)E(k) = E(k)$(m) for proportional 
angular rates. Finally, Appendices D and E contain computer programs 
for examples in Chapter IV. 



II. COMPUTATION OF THE FOUR PARAMETERS 

In this chapter, exact solutions for the four parameters when the 
angular rates are proportional to each other are derived. A numerical 
integration scheme is then selected for computational error analysis. 

CLOSED-FORM SOLUTION 

As shown in Chapter I, the time rate of change of the four para- 
meters is 

4(t) = ift(t)e(t) (II-l) 

“ 2 ~ 

where ei(t) is a 4 x 1 column matrix consisting of the four parameters 
e^, e- 2 * e 3 > e 4 an ^ ft(t) is a 4 x 4 skew-symmetric matrix of body angular 
rates as measured by the system gyroscopes 


0 -*z ~*y ~*x 

K 0 -* x K 

a(t) = . / (II-2) 

*y *X 0 -+z 

-*Y K 0 • 


A closed-form solution to (II-l) can be obtained if the angular rates 
are proportional to each other. Then the angular rate matrix J2(t) may 
be written in the following form: 

8(t) = Kf(t) (II-3) 
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where K is a constant 4x4 coefficient matrix defined by 


and f(t) is a scalar function. Under the above-mentioned assumption, 
the angular rate matrix at t^ and t 2 can be written as 


QCt x ) = Kf(t x ) 


0(t 2 ) = Kf (t 2 ) 


and is therefore commutative for all t. 


OC^) n(t 2 ) = a(t 2 ) fl(t ) 


The solution to (II-1) is given by [7, 8] 


e(t) = e 


1/ ft(x)dx 
- L o 


e(t 0 ) 


f<T)dT 


e(t 0 ) 


Let a(t) - / f(x) dx, 
t o 


then e_(t) - z 




, [It sW tt M + M t ...] .(t ) 

2 2 2 * 2 1 2 3 -3! ° 


( 11 - 10 ) 



8 


Since K is a skew-symmetric matrix, the following identities can be 
obtained* 

K 2 = -(k x 2 + k y 2 + k z 2 )l = - k 2 I 
K 3 = - k 2 K 
K 4 = k 4 I 


In general 


n-1 


K n = (-1) 2 k 11 " 1 K for n odd 


n 


= (-1) 2 k n I 


for n even 


( 11 - 11 ) 


where 


k 2 = k 2 +k 2 + k 2 
X y z 


( 11 - 12 ) 


Using these identities, equation (11-10) can be further simplified to 
e(t) = {I +-“^+ (^-) 2 (-k 2 I)(|,) + (5M.) 3 (_k 2 K)(|,) 


+ C^~-) 4 (k 4 DC|,) + (^-) 5 (k 4 K)(y,)...}e(t ) 


- ^ - C^-k) 2 1 + (SM k) 4 i - ...] 


+ |ja(|)k_ + _ ...J>e(t o ) 


= {l{Cos(^^)3 lSin(^|^)]>e(t o ) 


( 11 - 13 ) 
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Example 


Let - ajt 
*y ~ a 2 C 

K ■ V 


then 


K = 


k 2 = 


-a. 


a l + a 2 + a 3 


Sir 


-a. 


-a. 


and 


t (t) -‘/NdT- j- l^f 

o 2 o 2 


therefore 


2 2 

e^(t) = {I Cos )'+ ^ Sin(r|- )} e.(o) 


The same result is obtained by using the Peano-Baker method of successive 
approximation which is presented in Appendix A. 

Both equations (11-10) and (11-13) are exact solutions for the four 
parameters under the assumption that the angular rates are proportional 
to each other* Equation (11-13) is a closed-form solution which is 
obtained by making use of the fact that K is a skew-symmetric matrix. 
These exact solutions are expressed in terms of the angular rates* If 



10 


rate-integrating gyroscopes are used for the inertial system, then the 
outputs of the gyroscopes are the integrals of the input rates, i.e. 

ft " 

A0^ =J <f>. dt for i = x, y, z. 

fc o 

For this reason it will be necessary to express the exact solution in 
terms of the integral of the input angular rates. This can be developed 
in the following maimer. 

For proportional angular rates, the angular rotations about each 
coordiante axis can be represented by 

K. - k x f<t) 

* y - y <*> 

i z = kf(t) . (11-14) 

Therefore, the integral of the angular rates may be written as 
t 

A0^ = k ./ f(f)d(x) = k.a(t) for i = x, y and z. (11-15) 

x 1 c o 1 

Now expressing the arguments of equations (11-10) and (11-13) in terms 

of A0 , A0„ and A0 , the following expressions can be obtained, 
x y z 

a(t)k a(t)(k 2 + k 2 + k 2 )2 
_ x 2 z 

2 2 

2 2 2 xr 

(A0 X + A0y + A0 Z )^ 


2 


( 11 - 16 ) 
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K = a(t) 
k a(t)k 


0 

k 


-k 


-k 


-k 


-k 


-k 





k -k 

k 

0 





L x y 

z 


J 



0 

-A0„ 

-A0 

-A0 





z 

y 

X 



A0 

0 

-A0 

A0 


_ 1_ 


z 


X 

y- 

A0 


A0„ 

A3„ 

0 

-A0 




y 

X 


z 



A6 v 

-A0 

A0 

0 


_ A0 


X 

y 

z 


- 

A0 







and 







a(t)K . 

_ A0 





2 


2 





where (A0) 

2 

2 2 2 
= A0 + A0 + A0 






X y Z 




and 

r— 







0 

-A0 

-A0 

-A0 





z 

y 

X 



A3 

o • 

-A0 

A0 


A0 = 


z 


X 

y 



A0 

A0 

0 

-A0 




y 

X 


z 



A0 

-A0 

A0 

0 



L 

X 

y 

z 



Substitution of (11-16) , 

(11-17) 

and (11-18) 

yields 









A0 






(T"> 




e/t) 

=5 

e 

e(t ) 
— o 





and 


e(t) «' {I Cos A + M ' A6 


A6 Sin( “ )} £<*„> 


(11-17) 


(11-18) 


(11-13) 


(11-19) 


( 11 - 20 ) 
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Both equations (11-19) and (11-20) are exact solutions for the 
four parameters and are expressed in terms of the integrals of the input 
rates . 


Numerical -Integration Scheme 

As shoxm in Chapter I, the vector differential equation of the 
four parameters in terms of the body rates relative to the reference 
system is 

e.(t) = -j JKt)e.(t) 

Then e(t) = e(o) + — ft(t)e(t)dt (11-21). 

~ ~ o2 - 

where f z i fi(t)e(t)dt can be solved by various numerical integration 
0 2 “ 

techniques using a digital computer. A large number of numerical 
integration schemes have been proposed for the integration of this class 
of differential equations. The most commonly used integration schemes 
are the Euler algorithm, [4] 
where 

e[(n+l)T] = e[nT] + T e[nT] (11-22) 


the Modified Euler algorithm. 
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where 


e[(n+l)T] = e[nT] + £ {e[nT] + BIhLL (eJnT] + T e.[nT])} (11-23) 

“ z 2 

and the Fourth Order Range-Kutta algorithm, [3] 
where 

eJ(n+l)T] = e_[nT] + ~ {m^ + 2 ^ + 2m^ + m^} (11-24) 


=1 


= T e 


m = — . ft[nT]{e[nT] + im } 
“ 2 2 - 2 “1 

m„ = S2[nT]{e[nT] + — m_} 

—3 2 — 2 ^ 

m = y fl[nT].{e[nT] + m } 

• 2 3 


A different numerical integration scheme is considered in this study. 
This can be derived in the following manner. 

From equation (11-19), the exact solution for the four parameters is 


i(t ) - 4^ e(t ) 

— o 


= $(t,t )e(t ) 

O — Q 


(11-25) 
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where the matrix §(t,t o ) is called the state transition matrix. 


A0(t 2 ,t 1 ) + A0(t 1} t o ) A0(t 2 ,t 1 ) A0(t ls t o ) 

( ^ , o 2 

Since e ^ ^ = e e 


then 


< ^^ t 2’ t 0 ) $ ^ t 2’ t l^ ) ^l^cP for a11 fc 2» *]_» t 0 

This is the group property of the state transition matrix. From 

this group property, it is evident that for t = mT, T > 0, the recursive 

formula for equation (11-25) is 

, A0[(mKL)T, mTh 

^[(rntl)T] = s'- 2 ' eJmT] 

= $1 (m+l)T,mT]£[mT] , m = 0, 1,*«* (11-26) 

There are several different methods of evaluating the state 
transition matrix* The principal methods are [7]: 

1. The infinite series method 

2. The inverse Laplace transformation method 
3* The transfer function method 

4. The Sylvester’s theorem, and 

5. The Cayley-Hamilton technique 
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Among these methods, the infinite series method is most suitable for 

digital computation [16], 

(~) 

In the infinite series method, the state transition matrix e 
is calculated by the infinite series 


e 



!+(#■) + 


( f )2 + 
2 ! 


(M)3 

2 

3! 


(11-27) 


Since the infinite series (11-27) is uniformly convergent for all 
finite elements of A0[17], it can be computed by the truncated series 


„ P i 

$ = l 12/ (11-28) 

i=o * ! 


where A° = I 

within prescribed accuracy using a digital computer. Thus, for 
proportional angular rates, the numerical integration scheme for the 
vector differential equation of the four parameters is 

e[(mfl)T] = $ [ (mfl)T, mT]e[mT] , m = 0,1,.... (11-29) 

e(o) = e(o) 


where a hat O over a quantity denotes that quantity is an approximation 
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as a result of the finite series approximation. 

It has been shown by Marshall [18] that the Euler , the Modified 
Euler and the Fourth Order Range-Kutta algorithms are equivalent to 

A0 

the first two, three and five terms in the series expansion for e 2 ' , 

respectively. The authors also concluded that: 

(1) taking k terms in the series expansion of the matrix exponential 

series is equivalent to using a (k-1) st-order Range-Kutta numerical 
integration scheme. 

(2) a computer program written using the first k terms of the matrix 
exponential series will provide greater computational efficiency than a 
program written using a (k-1) st-order Range— Kutta numerical integration 
scheme. Therefore, by investigating the infinite series method, a large 
class of numerical integration schemes are being studied. 



III. COMPUTATION ERROR BOUNDS 

In Chapter II, a numerical integration scheme for solving the 
four parameter vector differential equation is derived. The numerical 
integration scheme will produce, corresponding to each mT, a vector 
£(mT) , which is an approximation to e^(mT) , the exact solution of the 
four parameters vector differential equation. The difference between 
e(mT) and <e (mT) is called the truncation error c(inT) * The truncation 
error is caused by the fact that only a finite number of terms of the 
infinite series is used in the numerical integration scheme. Due to 
the fact that all digital computers work with only a finite number of 
digits, the computed solution e*(mT) will in general not agree with 

A A A 

e_(mT) . The difference between e*(mT) and _e(mT) is called the roundoff 
error i:(mT) . This chapter analyzes both the truncation error and the 
roundoff error introduced in the floating-point computation of the 
four parameters using the finite series approximation method. Vector 
norms and matrix norms will be used to give an assessment of the size 
of a vector or a matrix, respectively. Their properties and definitions 
are given in Appendix B. 


Truncation Error 

As developed in Chapter II, the exact recursive formula for the 
four parameters is 


17 
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( AQ[(mfl)T] ) 

eJ(m+l)T] = e ^ _e[mT] 

= $ [ (m+1) T , mT]_e[mT] 


(III-l) 


for t = mT, T > 0 and m = 0,1, 2*»** where the matrix exponential 
( A0[(mfl)T] ) 

e'' ^ is defined by 

Ae[(nHhl)T3 oo ,AQ[(m+l)T],i - 

e 1 - l { 2 ) (III-2) 


and .AQ[ (mH-l)T] q_ 
v 2 y — X 


(IH-3) 


For digital computation, equation (III-l) is generated by the following 
approximate recursive formula: 


e[(m4-l)T] 


P 

l 

i=0 


,AG[ (m+1) ] - 

v 2 ' e [mT] 


i! 


= i[(nrt-l)T]e[mT] 


(III-4) 


The error incurred by using the approximate recursive formula will be 
considered for both constant angular rates and time varying angular 


rates . 
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Constant Angular Rates 

,A0[(mfl)T] v 

( 2 “) 

For constant angular rates, the matrix exponential e 

A0 

__ 

is a constant matrix e . The state transition matrix *[(irrKL)T, mT] 
is also a constant matrix which can be represented by 

A0 

$[(m+l)T, mT] = e ^ 



= $ for m = 0,1, • • • • (1II-5) 

where k is given by equation (II-4) . 

From equation (III-l) , the exact recursive formula for the four 
parameters is 

e[(uri*l)T] = 4>e[mT] (III-6) 

By a process of iteration, the following equation is obtained 

e[mT] = e(0) for m « 0,1, (III-7) 

where e^(O) is the initial condition of the four parameter vector and 
is defined to be the identity matrix. 
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Let $ = $ + E 


(III-8) 
<^> 


^ \ r\ 

where $ is the approximating matrix for the matrix exponential e ^ 


P Ar ) 1 

* = i c-4?-) 

i=0 i. 


and where E is the difference between the matrix exponential e 
and the approximating matrix $ 


(III-9) 




E = l 


K v i 

(oT) 


(III-10) 


i=p+l 1 • 


Substituting equation (III-8) into equation (III-7) yields 


e[mT] = [$ + E] m e(0) 


(III-ll) 


Expanding equation (III-ll) yields 

~rn ^rn— 1 — -f 

_e[mT] = 4> e^(O) + ( £ $ E $ •) e_(0) 

i=0 

+ 0(i m " 2 , E 2 ) e(0) (III-12) 

where 0(o m-2 , E 2 ) represents terms which are of higher order in E. 

If the four parameters are computed using the approximate recursive 


formula 
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e [ (ntKL)T] = * e[mT], (III-13) 

then the solution e(mT) to equation (III-13) is 

a A ni 

e[mT] = $ e(0) (III-14) 

The truncation error is defined as the difference between the 
exact solution e^[mT] and the approximate solution e[mT], Subtracting 
equation (111-14) from (111-12)' gives the propagation of the truncation 
error c (mT) 

c(mT) = e(mT) - e^(mT) 
m-1 

= ( l S 01 ” 1-1 E $ i )£(0) 
i=0 

+ 0($ Itl_2 , E 2 )e03) 

The norm of the truncation error is 

m-1 

I |c(mT> | | <_ || l i® -1 " 1 E * x | | * | |e(0) j ] 
i=0 

+ I |0(S m “ 2 , E 2 ) | j • ! |e(0) | | (III-15) 
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For practical purposes, the higher order terras in E are much 
smaller than the first term and therefore may be neglected yielding 


| c_(mT) | | < m * ||®!| ml * | |E| | * | |e(0) | | 


(III-16) 


f or l -L K l.i - Z i es s than one, the norm of E is given by [19] 

2(p+2) 


|| E || < llKl|P +1 : (t/ 2) p+1 . ( P +2) 

(p+1).' (p+2) - [ | K 1 | * (T/2)' 


(III-17) 


and it can be shown that 


nin"" 1 - ii ! ir 1 

i=0 

< e (m-l).||K||. T/2 _ (III-18) 


Substituting equations (III— 17) and (III—18) into equation (XIX— X6) gives 


||c(mT)|| < m • e (™-l)-llK||-T/2 . 


(III-19) 


(IMIZH : CT/2)P +1 
^ (P+1) ! 


(p+2) 

* (p+2) - I |Kj [ * (T/2) 


) • | | e(0) | | 


Equation (III-19) gives the truncation error bound in computing 
the four parameters using the finite series approximation method for 
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constant angular rates. .It is shown that for a fixed time period 
t = mT, the truncation error bound decreases with decreasing time 
increment T. For a fixed time increment T, the truncation error 
bound decreases with increasing number of terms used in the numerical 
approximation of the state transition matrix e KT /2 # 

Time Varying Angular Rates 

From equation (III-l) , the exact recursive formula for the four 
/ 

parameters is 

/ 

e[(m+l)T] = $[(nrt-l)T, mT]e[mT] (III-20) 

for t = raT, T > 0 and m = 0,1,2**** and from equation (III-13), the 
approximate recursive formula for the four parameters is 


e[(m+l)T] = $[ (nrhl)T, mT]e[mT] (III-21) 

for t = mT, T > 0 and m = 0,1,2* ••• 

By definition, the truncation error c^mT) is the difference between 
the exact solution and the approximate solution. Hence 

c_[(mhl)T] = e[(nri-l)T] - £[(nrKL)T] 


( 111 - 22 ) 
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To simplify the notation, let 

$>[(m+l)T, mT] = $(mfl) (III-23) 

$[(m+l)T, mT] = $(m+l) (III-24) 

and T be ommitted in [mT] 

Define the remainder matrix E (n+1) by 

E(mfl) = $(nri-l) - $(mfl) (III-25) 

V 

Now, substituting equations (III-20) and (III-21) into equation (III-22) 
yields 


c^m+1) = (nrt-1) £(m) - $(m+l)e^(m) (111-26) 


Utilizing equations (III-25) and (III-22) , equation (III-26) may be 

A A 

expressed in terms of $ (m+1) ,c (m) , E(m+1) and £(m) . Thus 


c(mfl) = $ (m+1) j;(m) + E(m+l)e^(m) + E (m+1) c^(m) 

for m « 0,1,2* * * * (III-27) 


Note that c(0) = ji(0) - e^(O) 
= 0 
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If the initial conditions £(0) are known, then equation (III-27) can 
be solved recursively for m = 0,1,2* ••• 

Thus, 


c(l) = E(l)e(0) 

£(2) = 8(2)E(l)e(0) + E(2)e(l) + E(2)E(1)£(0) 

c(3) = $(3)*(2)E(l)e(0) + $(3)E(2)e(l) + $(3)E(2)E(l)e(0) 

+ E(3)e(2) + E(3)$(2)E(1)£(0) + E(3)E(2)e(l) 

+ E(3)E(2)E('l)e(0) 

= $(3)S(2)E(l)e(0) + J(3)E(2)e(l) + E(3)e(2) 

+ o 3 (e 2 ) 


m-1 m- (i+2) ^ 

c^(m) = £ [ 7T $(m-k)] E(x+l)e^(i) 

1=0 k=0 


+ 0 m (E^ s i) for m = 1,2* (IXI-28) 


where O m (E 2 ,i) represents terms which are of higher order in E(m) and 
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ir <j>(m-k) = $(m)$(m-l) • • • • $(m-i+l) , i > 0 

k=0 

El , i = 0 (III-29) 


where I is the identity matrix. 

It follows from equation IU-2'8) that 

m-1 m-(i-2) 

I I c.( m ) | | < I [ it | | $ (m-k) | ] ] • |jE(i+l)|j • j |e(i) | | 

i=0 k=0 

+ l|0 m (E 2 ,$)|| (III-30) 

By using an approach similar to that used in the constant angular 
rate case in determining | | $ (m-k) | | and | |'E(i+l) j | , and by neglecting 
the higher order remainder terms a closed-form solution for the 
propagation of the truncation error bound can be obtained. Note that 
for constant angular rotations, equation (III-30) reduces to 

m-1 m-(i-2) 

l|c(m)[| < I [ w ||3j|] • | | E | | • ||e(i)|| 

i=0 k=0 

<m • I]i| r- 1 . I ] E I I • ||e(0)|| (III-31) 

This is in agreement with equation (III-16) which is derived by 
assuming constant angular rates. Observe that the higher order 

n A 

remainder terms 0 m (E ,$) in equation (III-28) are generated by the 
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product term E(m+l)e^(m) in'equation (III-27). Since for practical 
purposes, 

m-1 m-(i-2) 

[ ir | 1 $ (m-k) | | ] • l|E(i+l)|j . | |e(i) | | » || 
i=0 k=0 

the higher order remainder terms may be neglected. Consequently 
equation (HI-27) may be approximated by 

A ^ 
c^(nrKL) 53 $ (nrfl)_c (m) 4- E (m+1) e_(m) 

for m = 0,1,2 (111-32) 

Equation (111-32) gives the propagated truncation error which can be 
evaluated by means of a digital computer* 

An interesting form of solution for the propagated trunc ation error 
can be obtained by expressing equation (III-26) in terms of $(nrt-l), ciCm) , 
E(nH-l) and e(m)‘ 4 Thus, by utilizing equations (III-25) and (III-22) , 
equation (III-26) may be rewritten as 

c.(nrhL) = <£>(m+l) _c(m) + E(m4"l)£(m) - E(mH-l)c_(m) 

for m = 0,1,2' ••• (III-33) 

If the initial conditions e(Q) are known, then similar to equation 
(III-27), the solution to equation (III-33) is 


o (E 2 ,i) 1 1, 
m 
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m-1 m-(i+2) 

£(m) = l [ ir $(m-k)] E(i+l)e(i) - 0 (E ,$) 
i=0 k=0 


for m = 1,2* * • • 


(III-34) 


xjhere 0 m (E‘ i 


,$) represents terms which are higher order in E(m) and 


i-1 

ir O(m-k) = $ (m)$ (m-1) • ♦ ♦ • $(m-i+l) , i > 0 

k=0 

=1 , i = 0 (III-35) 


Using the fact that 


(1) e(i) - $(1)$ (i-1) • • • • $(l)e(0) 


i-1 

or je(i) = ir $(i-k)] e^(O) (III-36) 

k=0 


and (as shown in Appendix C) 


(2) for proportional angular rates 


0(m)E(k) - E(k)$(m) for all positive integers k and m 


equation (III-34) becomes 
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■m-1 m-(l+2) i-1 

c_(m) = £ E(i+1) [ it $(m-k)][ ir $(i-k)] e^(0) 

i=0 k=0 k=0 


- O m [E 2 ,#] 


for m = 1,2' 


m-1 m-1 

£(m) = £ E (i+1) [ tt $(m-k)] e_(0)' 

1=0 k=0 

k^(m-3.-l) 


- 0 [E^,*] 
m 


for m = 1,2* 


It follows from equation (III-37) that 


m-1 m-1 

_c( m ) | | .S l | | E(i+1) | | • [ ir | | $ (m-k) | | ] • | | e(0) | j 
i=0 k=0 

, k^ (m-i-1) 


+ I |0 m [E >$] | j 


for m = 1,'2* 


| ■ AOCx+1) ■ . 

For II 2 II less than one, the norm of E(i+1) satisfies 


||ECi + l>||< ||^|±iI||P+l 


(P+1) •' 


(P+2) 

, , A0(i+2) 
(p+2) - | 1 '2 ' 
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Equation (III-37) gives the propagated truncation error vector 
c^m) in computing the four parameters using the finite series 
approximation method for time-varying, proportional angular rates. 
Each element of _c (m) can be determined by neglecting the higher order 
remainder terms. Equation (III-38) gives the propagated truncation 
error norm. It shows that the truncation errors depend upon such 
factors as the initial conditions of the four parameters, the 
magnitude of the angular rotations and the number of terms, p+1, used 
in the numerical approximation of the state transition matrix. 
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Roundoff Error 

Generally there are two different approaches in analyzing the 
roundoff error in digital computation. These are the deterministic 
approach and the statistical approach. The deterministic approach is 
exemplified by Wilkerson’s work [20-26] on determining maximum bounds 
for the roundoff error. The statistical approach is advanced by 
Henriei [27-31] and has been verified by en tensive numerical 
experimentation. Since the truncation error bound derived in the last 
section is based on the deterministic approach* Wilkerson’s approach 
will be used in analyzing the roundoff error. The layout of this 
section is as follows. The roundoff error for the fundamental 
arithmetic operations will first be developed. Then the roundoff 
error bound in the computation of the four parameters using the finite 
series approximation method will be discussed. 

Roundoff Errors in Floating-Point Computation [20*32] 

111.1 Notation 

For any real number x, let x* be its floating-point machine 
representation. For floating-point computations* let fl[.] be the 
floating-point machine number obtained by performing the arithmetic 
operation specified by the parenthesis [7]. It is assumed that 
computation proceeds from left to right. 

111. 2 Floating-Point Machine Number Representation 

In floating-point arithmetic, all numbers are represented in the 
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computer by floating-point machine numbers which are of the form: 

x* = (sign x) * • a ' (III-40) 

where a is a terminating g-nary fraction satisfying the following 
normalization condition 

1 <_ a < 1 (III-41) 

g 

b is an integer, ranging between -E to E, and g is the base of the 
number system employed by the computer. The g-nary fraction a is called 
mantissa or the fractional part of the floating-point machine number x* 

It is represented by 

t 

a = l a.g -1 (III-42) 

1-1 

where t is the number of g-nary digits a computer used for the 
fractional part. The integer b is called the’ exponent or the charac- 
ter which is given by: 

b = [logp[x| ] + 1 (III-43) 

where the brackets [♦] denote the largest integer not exceeding the 
quantity inside the brackets, and log^ denotes the logarithm to the base g. 
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The range covered by the magnitude of a floating-point machine 
number x* is 

(fT 1 )!? -33 < |x*| < (1 - 3 -t ) • g E 

or g"(E+D < | x *j < (1 - 3 -t ) • g E (III-44) 

The range of the computer is defined by the interval 

R = [-(1 - $ _t ) ■ f$ E , (1 - p“ c ) • g E ] (XII-45) 

It is assumed that enough bits are allowed for the exponent so that 
no computed floating-point machine number will lie outside the 
permissible range ♦ 

III. 3 Input' Roundoff Error 

Consider a real number x. The process of replacing x by a floating- 
point machine number x* is called input rounding. Input rounding can 
usually be achieved either by truncation or by rounding. In truncation, 
the first t digits of the mantissa are retained and those digits 

beyond the first t digits are dropped. Since x* = 

r 

.(sign x)'3 b *( £ a.3 -i ), it is evident that if x e R, |x| 
i=l 1 

the input roundoff error is bounded by 
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jx - x*| <_ (3 b • 6 -t (IIX-46) 

Noting that 3 ^ | x | ♦ equation (III-46) becomes 

|x - x*| < |x| • 3 1-t (HI-47) 

or fl[x] ® x* 

» x(l + e, ) (III-48) 

in 

where |e in | £ (III-49) 

The above relation shows that if x e R and | x | g (E+l) ^ t | ien 
the relative error of the truncated floating-point representation x* 
of x is at most 

In roundings the first t digits of the mantissa are retained 
after a 3/2 is added to the (t + l)th digit. Therefore the input 
roundoff error is bounded by 

|x - x*| < 6 b • (| • r (t+1) ). (IIX-50) 

Since 3^ <_ |xj ’ 3, equation (XII-50) may be expressed as 


jx - x*| <_ |x| • ~ • 3 1 


(III-51) 
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or fl[xl = x* 

« x(l + e in ) (III-52). 

where Is. | < i • (III-53) 

1 in 1 — 2 

The above relation shows that if x e R and |x| >_ 6 (E+l)» t ^ en 
the relative error of the rounded floating-point representation x* 
of x is at most ^ 

III. 4 Addition and Subtraction 

Consider the addition or subtraction of two floating-point machine 
numbers x* and y* each with a t digit mantissa. Let 

b t 

x* = (sign x) ’ ft x • ( £ a . S 1 ) (III-54) 

i=l X ’ 

b t 

y* = (sign y) • 6 y • ( £ a .g -1 ) (III-55) 

i=l y ’ 


b 


X 


< 


b 


y 


b 

and fl[x* ± y*] = z* = (sign z) • B Z 


t 


a 

i=l 



(III-56) 


It is assumed that the sum (or difference) of x* and y* is computed 

in the following manner. The exponents b and b are compared and the 

x y 

fraction of y* is right-shifted b - b places. The fractions are 

A y 

then added algebraically to form an intermediate sum IS. This 
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intermediate sum consists .of. (t + 1) digits and a possible carry. 

The extra digit is a guard digit obtained from the fraction which is 
shifted right. After the addition, the intermediate sum is left 
shifted or right shifted so that the resulting mantissa satisfies the 
normalization condition, the exponent b x being adjusted accordingly. 
Finally the resulting mantissa is truncated or rounded to t digits. 
This gives a . Rounding by truncation will be assumed from here on. 

The process may be illustrated by three examples of addition of 
machine numbers in 4 digits floating-point decimal arithmetic. 

Example 1: — < IS < 1 

g - 

fltio 4 (0.7414) + 10 1 (0.3995)] = 10 4 (0.7417) 

ay is shifted 3 spaces to the right and the addition takes place 
in the form 

guard digit 

10 4 X 0.7414 0 
-(-10 4 X 0.0003 9 
10 4 X 0.7417 9 

The intermediate sum is 10 4 X 0.74179 which is then normalized and 
truncated to 10 4 X .7417. 

Example 2 : IS > 1 

fl[10 5 (0.7419) + 105(0.6159)] = 10® (0.1357) 

The addition takes place in the form 
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guard digit 

10 5 X 0.7419 0 
+10 5 X 0.6159 0 
10 5 X 1.3578 0 

The intermediate sum is 10^ X 1.35780 which is then normalized and 
truncated to 10^ X .1357. 

Example 3 : IS<-g- 

fl[10 _4 (.1000) + 10 -6 (-.9999) ] = 10" 5 (.9001) for truncation 
The addition takes place in the form 

guard digit 

10" 4 X .1000 0 
-10~ 4 X .0099 9 
10 -4 X .0900 1 


The intermediate sum is 10 4 X .09001 which is normalized and 
truncated to 10 - "* X .9001. 

b t 

If the computed sum is (sign z) • 3 • ( V a ,3 ) 

i=i ^ 


then it is evident that the magnitude of the error is bounded by 


b 

|(x* ± y*) - fl[ x * ± y*] | £ p Z ■ p 


(III-57) 
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-Z 

Since 8 <_ ] (x* ± y*) | • 8 for (x* ± y*) 4 0 , 

equation (III-57) may be rewritten as 

] (x* ± y*) — fllx* ± y*] J 1 (x* + y*) | • g^ 
or fl£x* ± y*J = (x* + y*) (1 + 6) (III-58) 


where 

| <S | g 1 " 1 (III-59) 


Thus the relative error of the truncated sum (or difference) of x* and 
y* is at most B^*' t . 

III. 5 Multiplication 

Consider the multiplication of two floating-point machine numbers 
x* and y* each with a t digit mantissa. Let 


b x t 

x* = (sign x) * 6 * ( I a 8 X ) 

1-1 


b t 

y* = (sign y) • 8 y • (- J a -g -1 ) 

i=l 


and fllx* X y*J - P* = (sign P) * B P ( I a . 8 -i ) 

i=l ^ ,:L 


(III-60) 

(III-61) 

(III-62) 


It is assumed that the product of x* and y* is computed in the 
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following manner. The exponents b^ and are added together and the 

product of t . t 

£ a . 6 1 and J a . 6 1 is then computed. The 
i=l 1=1 y,± 

resulting intermediate product will have a fractional part of 2t or 
(2t - 1) digits. This product is normalized if necessary by a left- 
shift, the exponent being adjusted accordingly. The resulting product 
is then truncated to give a t digit mantissa of the computed product P*. 

Example 

f 11.1303 X . 1003] = 10 -1 X .1306 

Absolute error = |C»1303 X .1003 - fl[.1303 X . 1003] ] = .909 X 10“ 5 < 10 -5 

Relative error = | (.1303 X .1003) - fll.1303 X .1003] | „ , rt _3 

(.1303 X .1003) * 696 X 10 

If P* e R, then it is evident that the magnitude of the roundoff 

error is bounded by 

b p 

I (x* X y*) - fljx* X y*] I < $ • e" E • (III-63) 

bp 

Since 3 |x* X y*| * 3> equation (III-63) may be expressed as 

| (x* X y*) - fllx* X y*] ] £ |x* X y*| . 3 1-t 
or fllx* X y*] = (x* X y*) (1 + g) (III-64) 


where | \ \ <_ 3 1_t 


(III-65) 
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Thus the relative error of the truncated product of x* and y* 
is at most 
III. 6 Division 

Consider the division of two floating-point machine numbers x* and 
y* each with a t digit mantissa. Let 


b t 

c* = Csign x) * 3 X ’ ( ^ a . B -i ) 

i=l X » 1 


by t 

y* = (sign y) * e • ( I a yj± B -1 ) ± 0 


and fllx* f y*]= D* = (sign D) • 3 * ( Z a . 

i=l D,;L 


(III-66) 

(III-67) 

(III-68) 
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Absolute error =| (10 - ^ X .9137 ^ 10 2 X .1312) 

- fl‘110" 6 X .9137 * 10“ 2 X . 1312] 1 

= .17 X 10“ 7 < 10~ 7 

Relative error = Absolute error 

(10 6 X .9137 * 10“ 2 X .1312) 

= .24 X 10” 4 < 10 -3 

If D* e R, that it is evident that the magnitude of the roundoff 
error for division is bounded by 

b 

] (x* * y*) - fllx* ~ y*] <. 8 ^ * B _t (III-69) 

bjj 

Since 3 <_ ]x* * y*| • 6 > equation (III-69) becomes 

j (x* 4- y*) - f ljx* 4 y * } | <_ ]x* ^ y*| ' / 

or fljx* v y*J = (x* •=* y*) (1 + q) ~ (III-70) 

where ]n] ^ 6 1 " t 

Thus the relative error of the truncated quotient of x* and y* 
is at most for y* not equal to zero, 

III . 7 Extended Additions 

Consider the addition of a sequence of n floating-point machine 
numbers x^*, x 9 *,**** 5 x* . 

Z. T1 
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Let Sj = fl&Jj 



and S* = £l£S*_^ + x*] for i > 1 (HI-71) 

Then by applying equation (III-58) to equation (III— 71) » the computed 
sum for the first two terms of the sequence can be represented as 

sl = flfx* + x *] = x*(l + fi 2 ) + x*(l + « 2 ) 

where |fi 2 | <_ g 1 "* (III-72) 

Similarly the computed sum for the first three terms of the sequence 
can be written as 

S 3 = fUS* + **] 

- S*( 1 + 63 ) + x*(l + 5 3 ) (III-73) 

where ] 6 ^] <_ g 1 C 

Substituting equation (III-72) into equation (III-73) yields 


S* = x*(l + 6 2 )(1 + 63 ) + x*(l + & 3 ) 


(III-74) 
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It follows that the computed sum for the sequence of n terms can be 
represented as 


n 


fl[s Li + 4 

x-^ Cl + 62) (1 + S3) ••*•(! +* 6 n ) + 

xjjCl + 63) (1 + + 6 n ) + 

x (1 + 6 ) (1 + 5, ) * • • • (1 + 6 n ) + * • » • + 

j j ^ 


x (1 + 6 ) 
n n 


(III-74) 


where ] 6 1 <_ for 1 - 2 , #### n«. 

Expression (III-74) shows that the upper bound for the roundoff 
error is least when the smallest terras are added first, since the 
largest factor, (1 + <$ 3 ) 0- + S 3 )* # ‘*(l 4- <$ n ) , is associated with the 
smallest term. 

III. 8 Extended Product 

Consider the multiplication of a sequence of n floating-point 
machine numbers x-^ a x 2 5 x 3 > * • # * > x n * 


Let p* = fl[x^] 


= x 
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and p* = flip? . .xlh for i > 1 (III-75) 

i ^ 1 

Then by applying equation (IIIr-64) to equation (JII*«75) , the computed 
product for the first two terms of the sequence can be represented as 

P 2 - fl[pi X^] 

= x? x*(l + ? 2 ) (III-76) 

where ]s 2 | 1. g 1_t 


Similarly the computed product for the sequence of n terms can be 
expressed as 


p ' = tp. 


n—l 


x ] 
n 


X 1 x 2 + + 


(1 + £ n ) 


where ] <_ g 


1-t 


for i — 2, • * * *n. 


(III-77) 


The actual error incurred will depend on the order in which the 
multiplications are computed, but the error bound given by equation 
(III-77) is independent of the order of multiplication. 

III. 9 Roundoff Error in Matrix Operations 

Based upon the previous derived error bounds, it can be shown that 
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if k is a scalar, A and B are n X n matrices, then 


'* ‘ + Sin, «>J 


(III-78) 


flU* + B*] = Ka^ + + 6 )] 


(III-79) 


fllk* A*] = Jk*a. . (1 +4..)] 
H xj 


(III-80) 


A * 

where a . , a . . and b „ . denote the (i , i ) element of the matrices 
ij ij ij ,J 

A, A* and B respectively. The e. T s are in general different but 

ij 

all are bounded by g t . The same is true for the S. *s and 

ij 3 .J 

For matrix multiplication, consider the multiplication of two 
n X n matrices A* and B* with elements that are floating-point machine 
numbers. Let c*. be the (i,j) element of A*B* which can be represented 
by 


fllcjj = flla^ b*. + a * 2 b*. + 


b*. ] 


xn nj 


(III-81) 


By applying equations (IXI-64) and (III-74) , equation (XII-81) becomes 

fllc*j] = la^b* (1+C 1 )(1 + « 2 )....(1 + 6 n ) + 

a i2 b 2j (1 + ? 2 ) (1 + V * * “ (1 + 6 n } + 
a*^b* . {1 + X > (1 + 63) •••*(! + 63) + .... + 
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a b . (1 + ? )(! + <$ )j 



in nj 

n 

n 


where 

itiii 

g i- t 


i = 

and 

i«ii i 

gi-f 


i = 


(III-82) 

■Lj « • * »Tl 

2, . . . .n 


Equation (III-82) can be written as 


- 1(1 + «„) ^ a* k bj.] 

n 

for £ a* k b^ ^ 0 and fl{c*J ^ 0 
k-1 

where 

a-6 ) n <_ 1 + a 1 (1 + 8 1-t ) n 


(III-83) 


(III-84) 


Since the roundoff error bounds to be derived do not depend 
critically upon whether equation (III-82) or equation (III-83) is used, 
equation (III-83) is assumed without loss in generality* The error 
bound for the last expression is very conservative, but it greatly 
simplifies the derivation of roundoff error bounds for extended 
matrix operation* 

Wow consider the roundoff error made in raising a n X n matrix A* 
to its p^h- poxtfer where p is a positive integer. Consider first the 
computation of A*A*. Let a” be the (i,j) element of A*2. From 
equation (III-83) , the computed value of 

ij 


can be represented as 
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fu.;®] 


I Cl + a (1) ) a*( 2 h 
ij 3.J 


(III-85) 


where 


(i - g 1-t ) < (1 + af-P) 11 < (1 + $ 1 ~ t ) n 

~ ij “ 


(III-86) 


Consider next the computation of A*fllA*2] . Let a* ^ ^ be the 

ij 

(i,j) element of A x 3 . From equation (III-82) , the computed value of 


a? can be represented as 


fHa*P>] = la a*f 2) (l + o (1) )(l + 0(1 + 6 ) 

13 ll lj 1 ^ 1 ^ 


a* a*f 2 > (1 + o VJ -')(l + OU +6) (1 + 5 ) + 

i2 2] 2j 2 2 n 

a* a* (D (1 + a^)(l + X, ) (1 + 8,) V 


lj 

.( 1 ) 


(1 + 5 ) + 


n 


13 


3j 


3j 

Cl), 


a* a*< 2) (l + a v 7)Cl + 0(1 + 6 )J 
xn nj nj 3 n 


(1+5 ) + 

n 


(III -87) 


where (1 - g'*''' t ) n 1 + (1 + g"*" t ) n 


]C ± | <_^ 
\* ± \ 


x L j • • • * xx 1 
x = l^*** *n « 


(III-88) 


Tor the computation of roundoff error bounds, equation (III-87) can be 


rewritten as 
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flIa AC3) J = 10 . + 


a 0) )a *»)j 

“ij )a ±3 J 


where (1 ^ 1 + cif. . ^(14- fi 1 "'**) 

j-j - 

Similarly , . it can Be shown that if p is an integer 


(III-89) 

\ 


flla*^ - ia + 1) )a* (p) J 

ij il ij 


(III— 90) 


. n ol-tx Cp-l) n ' (p-1) ' _l-t. Cp-l)n 

where (1—6 ) 1 + <_ (1 + B ) 


&II-91) 


Roundoff Error in the Computation of "the four Parameters 

Now consider the roundoff error incurred in the floatingpoint 
computation of the four parameters using the approximate recursive 
formula. From equation (III— 4) and using the same rotation as 
defined by equation (III-24) , the theoretical approximate recursive 
formula is 


p(.Ae J (m+l)] si 

£(m+l) = ( 1 -77 ) £(m) 

i=0 *• 


= S(m+l)e(m) (III-92) 

To compute e_(nrt-l) } the theoretical approximations 8 (nri-1) and e/m) 
are computed first giving the computed approximations § (nri-1) and e*(m) 
respectively. Then 8*(m+l) is multiplied by e* (m) to give £ (nH-l) 
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Hence the computed values of js (m+1) can be represented as 

e*(mKL) - fl[$*(ntKL)e*(m) ] (III-93) 

Let <j> . , (nri-1) be the (i,j) element of (nrKL) and e?(m) be the ith 
element of (m). The computed value of e^fritf-l) ,e* (m-KL) , can be 
represented as 

e*(mfl) ® ^(m+lJe^mHl + C ) (1 + 6 2 >(1 + 63 ) (1 + 64 ) + 

^(mrHOe^m) (l + 5^) (1 + <$ 2 ) (1 + 63) (1 + 64) + 

*i3^ +1 ^3^ m ) ^ + ^3) 63) (1 + fi ) 4* 

& 

^(m+DeJWd + (1 + § ) (III-94) 

Hence 

e^(m+l) = (ntH) e~ (m) (1 + o^) 

4>* 2 <Brf-l)e*(m) (i + 0 . 2 ) 

4 * (mfl)e*Cm)(l + ff .J 

1 3 d 10 

<j> (m+l)e (m)(l + a ) 

14 4 1 4 


(III-95) 
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where (1 — <^,1 + cr_^ <__ (1 + £-*- t ) <!> 

(1 - a 1 " 41 ) 6 '^ < 1 + a. . < (1 + g 1 ^) 6 ^ (j = 2,3,4) (III-96) 


Therefore, associating the factor (1 + o„) with. <f> ^ (m+1) , equation 
(III— 93) can be rewritten as 


e*(m+l) = [$*.(m+l) • (1 + a,,)] e*(m) 
— ij 13 - 

It will be shown in section 111*10 that 


(III-97) 


$*.(m+l) • (1 + a..) = <j>. .(m+1) + r. .(m+1) (111-98) 

13 13 13 13 

where r^ (nri-1) is called the local roundoff error. Substituting 
equation (III-98) into equation (III-97) yields 

e* (mi-1) = [$. . (m+L) + r. . (m+1) ]e* (m) (III-99) 

13 i 3 

Let R(m+1) = Ir^j (m+1) ] , then equation (III-99) can be rewritten as 

e* (mfl) = l8 (m+1) + R(m+l)]e* (m) (III-100) 

By a process of iterations, the solution On) for equation ( 111 - 100 ) 
is given by 
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^ lll J. . 

e (m) = ir [§(m-k) + R(m-k)] e* (0) (III-101) 

lc=0 

where 

nhJ. „ 

tt lS(jn-4c) + R(m-k)J = Ii’(m) + R(m) 1) + R(m-1)].... 
k=0 

... I$(l) + R(l)3 , m > 0 

sl , to = 0 (III-102) 

If the initial values of jj(0) are equal to the floating-point machine 
values, then equation (III— 101) can be written as 

_ to— 1 

e* (m) = { u li(m-k) + R(m-k)]}e(0) (III-103) 

k=0 “ 


Expanding equation (III-103) yields 


m-1^ m-1 m-(i+2) A i-1^ 

e* (m) ={ ir $(m-k) + J { ir $(m-k)] ♦ R(i+1) • I it 8(i-k)] 
k-0 i=0 k=0 k=0 


+ 0 m lR2$J } - e_(0) 


(III-104) 


Equation (III-104) gives the computed solution for the four 
parameters using the finite series approximation method. The theoretical 
approximate solution e^(m) for the four parameters can be determined 
from equation (III-92) by the same process of interation used to obtain 
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equation (III-101) . Thus 
m-sl^ 

e_(m) = { tt $(ra-^k) > •' e(0) (111-105) 

k=0 

The difference between the computed solution e* (m) . and the theoretical 
approximate solution e/m) is defined as the accumulated roundoff 
error . Hence, by subtracting equation (III— 105) from equation (111-104), 
the accumulated roundoff error :r(m) is obtained. Thus 

r/m) - e* (m) — e_(m) 

m^l m— (i+2)^ i— 

= £ X $(m-5k)j • R(i-HL) • l tt $(i-k)] • e/0) 

i=0 k=0 k=0 

+ 0 m lR?$J • e(0) (III-106) 

It follows from equation (III— 106) that the norm of the accumulated 
roundoff error is given by 

m— 1 m-1 

H^ Cm) li 1 J 0 ll R(i+1) !l * I^HiCm-k)]]] • ||e(0)|| 

k=^ (m-i-1) 

+ ||0 m [R 2 ,$]| j • ||e(0)|| for m = 1,2 (III-107) 


For practical purposes. 
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m-1 m-1 

l |jR(i+l)||[ tt | [ S (m-k) j | ] * |[e(0)|| » | |O m [R2,*]| | * | | e(0) j | 

i=0 k=0 

(m-i-1) 

Therefore, the higher order terms in R may be neglected yielding 
m-1 m-1 

| |r(®)| | 1 l | |R(i+l)| ! ‘ [ it |l$(m-k)||] • j|jeC 0 )|| 
i=0 k=0 

(m-i-0) 

for m = 1,2,.... (III-108) 

Equation (III-108) gives the accumulated roundoff error norm in 
computing the four parameters using the finite series method. It 
shows that the roundoff errors depend upon such factors as the 
initial conditions of the four parameters, the magnitude of the 
angular rotations and the local roundoff errors. 

Ill ‘10 An Example of the Procedures for Bounding the Roundoff Error 
Norm r(m) 

Consider the computation of the four parameters for constant 

angular rates. From equation (III-8) , the approximate state transition 

^ KT / 2 

matrix $ for the matrix exponential e is 

$ = I + KT/2 + (KT/2) 2 f 2 + (KT/2) p f p (III-109) 
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1 (j>\ 

where = — t for i ” 2,3,. ••• p« Let ^-±a j be the (i,j) element 
x « J 

of K^, Then the (i,j) element of 0 can be represented as 

L, = A ±i + k^T/2 + k^?)(T/2) 2 f, + + k^ (T/2)^f „ -§- 

xj xj z - L J £ 

+ k£p (T/2)Pf (III-110) 

where A-jj = 1 for i = j 

= 0 for i f j 

The application of equations (III-48) , (III-77) and (III-90) .leads to 

f£[k$p] = (1 + a )k£P 

ij 1 3 

f£[T] = T(1 + e. T ) (III-lll) 

xn j i 

fill 1 ] = T^(l + e £) 

T 

mi' 1 ] = 2^(1 + e _) 

"A 

and f.£[f £ ] = f £ (l + e ) 

f l 
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l_t 5£-4 (l- 1) 1-t 5£-4 

where (1. - g ) < 1 + a. . < (1 + g ) 

- ij - 


(i - g 1 ^) < i + e . < a + e 1 fc ) 

m — 


(III-112) 


n „ 1-tv 26-1 , , n , l-t.2l-l 

(1 - $ ) <_ 1 + <_ (1 + g ) 


(1 - g r ) < 1 + e f <(1 + g ) 

2T ~ 


and <1 - e 1 " 1 :)' < 1 + c, < (1 + 6 1_t > 

h~ 


Hence 


ark“<T/2)^] = k<f ( T /2)\ (i + , >t ) 


(III-13) 


where (1 - g 1_t ) 8 ^ 2 < 1 + e*. p < (1+ g 1 ^) 8 ^ -2 

— * IJ 


(iri-ii4) 


Now, the computed value of ^ can be determined by applying equations 
(III-74) and (III-113) to equation (III-110) . Thus 


♦ y - 4 y (l + kp) + k" 1/2 (1 + ^ + k^(T/2) 2 £ 2 (l + e p _ 2 )+ 


/m /o\ ^ / 


+ + vt> + — 
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+ (T/2) p f (1 + £ ) (HI-115) 

1J P 0 

where (1 - g 1 t ) p 1 + <_ (1 + g 1-t ) p 

(i _ < i + £p _ 1 £ (1 + g 1 ^)^ 4 

and (1 - g^ -t )^“^ + P <l + e -,<(1 + gl - t)7£-l+p f or £ - 2 , 3 ,....p. 

— p-x. — 

It follows that 

( 1 ) 

hi a + °ij ) = i ij (1 + S } + k ij T/2(1 + S-l } + *' • ' 

+ k[ 2 ) (T/2) 2 f 9 (l + 5 J + .... 

^ p -2 

+ k'f (T/2) \(1 + W + ‘ 

+ k//(T/2) p f (1 + c ) (III-116) 

u 

where (1 + g 1 “ t )P + ^"3 £ 1 + C p £ (1 + 1 3 1 -t )P +6 "j (III-117) 

(1 + g 1-t ) < ! + , < (1 + e 1 ’ 6 ) (P + «+ 6 -j (IH-H 8 ) 
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and (1 + g 1 "*) < 7 £- 1+ P) +6 -a < i + 


C . < (1 + 6 1_t ) ( 7 £~l+P)+6-j 

p ” / C 


j - 2,3,4. (III-119) 


Now, if i is an integer and i • g 1-t < • 1, then 

Cl - g 1- V < 1 + c < (1 + 


may be replaced by the simpler inequality [20] 


| c| < i O (III-120) 

where a = 1.06 g^ C 

Therefore, inequalities (III-117) , (III-118) , and (III-119) maybe 
replaced by the following inequalities 


.Spl 1 (p+6-6) a 


(III-121) 


ICp.il 1 (p+10-j) o 


(III-122) 


|c p -l\ 1 ( 7£+5+p- j ) cf 


(III-123) 
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Equation (III-116) may be rewritten as 

i?.(l + o..) = A ±i + kJP T/2 + kff (T/2) 2 f 0 + .... 
ij ij J -J X J ij 2 

+ k£p (T/2) £ f £ + . . . . 

+ k£p (T/2) p f p + r ±j 

/V 

= *ij + r ij (III-124) 

where the local roundoff error r_jj is bounded by 

l r ij! 1 a[(p+6-j)A ij (p+10-j)|k^ ) ! T/2 + (19+p-j) | k^ 2) j (T/2) 2 f 2 +. . . . 

ro\ 

. . . .+ (7£k5+p-j)| k | (T/2) £ f £ +. . . . 

+ (8p+5-j)|ki?*| (T/2) p f ] (III-125) 

J P 

Since j 2, inequality (III-125) may be rewritten as 

| r ij|<o[Cp+4)A ij + Cp+8)| k^j^J (T/2) + (p+17)| k£P| (T/2) 2 f 2 + 

.... + (p+7 +3)|k£P| (T/2)^f £ + .... 


.... + (8p+3)| k|p)| (T/2)Pfp] 


(111-126) 
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Let |R| denote the local roundoff error matrix with elements [ r^ ^ | , then 
|R| <_a[(p+4)I+ (p+8) | K j (T/2) + (p+17) | K | 2 (T/2) 2 f 2 + .... 

+ (p+7£+3) |K|£(T/2)^f + 

.... + (8p+3)|KjP(T/2)Pf p ] (III-127) 

It follows that 

| | Rj | £o[(p+4)l+ (p+ 8) [ j K| | (T/2) + (p+17)| |Kj | 2 (T/2) 2 f 2 + .... 
.... + (p+7£+3)| | K| |^(T/2)^f£+ 

.... + (8p+3)| |K| | p (T/2) p f p ] (III-128) 

Since £ | | k] | ^(T/2)^f . <_ e \ it can be shown that < 
i=0 

| | R| | < a [ (p+3) 1 I K 1 |T/2 _j_ i_ 2 j j K ] | t/ 2+7| j Kj ] (T/2) M K 1 I T/2-j 


(III-129) 
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Using the fact that for constant angular rate 


YnMlHm-WlI] <m t <»-U||K||l/2 

1=0 k=0 

(m-i-1) 


(III-130) 


and substituting inequalities (III-129) and (111-130) into inequality 
(III- 108) , the accumulated roundoff error norm j |r(m) | | bound is 


|r(m) | 


I K I l T//2 • a [ (p+3)e ^ i T ^ 2 + 1—2 1 |k| |t/ 2 
+ 7||K||Cr/2)J! K l l T/2 ]. | I e.(0) 1 | (III-131) 



IV. STUDY RESULTS 


• In order to check the validity and to demonstrate the applicability of 
the analytical results developed in the preceeding chapters, two exam- 
ples will be considered in this chapter. 

Example 1 

Consider the following first order linear fixed autonomous system 
x = Ax (IV-1) 


where x is a two dimensional column vector and A is a constant 2x2 
matrix given by 


A = 




(IV-2) 


Let the initial conditions for equation (IV-1) be specified as 


3^(0) 


1 

[ 

X 2 (0) 


1 

t— J 

- i 


(IV-3) 


The solution of equation (IV-1) at t = 1 sec. is to be computed using 
the following recursive formula; 
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x[mT] = e AT x! (m - 1)T] (lV-=4) 

for mT = 1, T > 0 and m = 1,2..., and e A ^ is to be computed by the 
truncated infinite series 

e A = l (IV— 5) 

i=0 11 


It is desired to determine the actual computational error norms and the 
theoretical computational error norms of x(l) so that the two error 
norms can be compared* 

Actual Computational Error Norm 

A digital-computer program is written to compute the actual com- 
putational error norms* The program is written in FORTRAN IV and has 
been run successfully on the IBM 360/50 digital computer at Auburn Com- 
puter Center, Auburn, Alabama. The actual computational errors are 
taken as the difference between the computed values of x(0 and the 
theoretical values of x(0 • The computed values of x(t) are obtained 
by implementing equation (IV-4) and equation (IV-5). The values of 
x(l) are then computed for p = 6,8 and T = 2^ ~ i = 1,2,... 10, 
using single precision (six hexadecimal digits or six bytes). The 
theoretical values of x(l) are determined by using the Laplace Trans- 
formation method. They are given by 
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x^l) 

x 2 Cl) 


(2e -1 - e~ 2 ) + (e" 1 - e 2 ) 
2 (e" 2 - e - " 1 ) + (2e~ 2 - e -1 ) 


and are computed in double precision (14 hexadecimal digits) * The re- 
sulting actual computational error norms are plotted in Pig* 1 as a 
function of the time increment T for p “ 6 and 8. Note the shape of 
the characteristic curve of the actual computation error ’norm* It is 
observed that minimum computational error norm occurs at T = *125 
second and T = .25 second for p = 6 and p = 8, respectively. 

Theoretical Computational Error Norm 

Prom equation tHI^19), the norm of the truncation error is bounded 
by 

(m-l)-||A||-T || A ||f + V + 1 _ p + 2 

(p + 1)! (p+2)! - I |a| |T 

• I k(°) 1 1 


Following a similar technique used in deriving equation (IXI-131) and 
noting that T, n! , and the elements of A are machine numbers, it can 
be shown that the norm of the roundoff error is bounded by 



Error Norm 



Figure 1. Actual computational error norms as a function of time 

increment for p=6 and p=8 on a 6-bytes fractional computer 
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m*e 


< "'- 1 >-|WI%JWI I + om JI A l! T - 1 >- ll A ll T J . | |A l|W 


+ 2ctIX + | |a| )Te^ A I |x(Oj I j 


The theoretical computational error norm is then the sum of the trunca- 
tion and roundoff error norm. The three norms are computed for p - 6 
and 8. The resulting error norms are plotted in Fig. 2 as a function 
of the time increment for p = 6 and 8. 

From Fig. 2, it may be seen that minimum theoretical computational 
error occurs at T = .125 second and T = .25 second for p = 6 and p - 8, 
respectively. This is in good agreement with the experimental results. 
Note also that for T greater than the optimal T, the computational error 
is dominated by the truncation error and the roundoff error can be ne- 
glected. For T less than the optimal T, the computational error is dom- 
inated by the roundoff error and the truncation error can be ignored. 
This shows that there are essentially two regions of computational er- 
ror. These are, due to their origin, the truncation region and the 
roundoff region. 

Fig. 2 also illustrates that, in the truncation region, the com- 
putation error is a function of both the time increment T and the order 
of the finite series p. Decreasing the time increment decreases the 
computational error. Increasing the order of the finite series p also 
reduces the computational error and increases the slope of the trunca- 
tion curve. In the roundoff region, the computation error is also a 
function of both time increment T and the order of the finite series p. 
Increasing the time increment results in a lower computational error. 




,01 a i 

T(sec) 


Figure 2. Theoretical truncation, roundoff and computational error norms 
as a function of time increment for p=6 and p=8. 
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Decreasing the order of the finite series also decreases the computa- 
tional error* Notice also that the slope of the roundoff error line is 
approximately minus one which compares favorably with the experimental 
results* 

The theoretical computational error curve is compared with the 
actual computational error curve in Pig* 3 for p - 8* It shows that the 
theoretical error norm is larger than the actual error norm. This will 
always be true since the theoretical result is an upper bound on the 
error. 


Example 2 

To check the theoretical results derived in Chapter III, the 
floating-point computation of the four parameters using the finite se- 
ries method is considered. 

Actual Computational Error 

A digital- computer program is written to compute the actual com- 
putational error* The computed values of e_(t) are obtained by imple- 
menting equation (11-29) for angular rates of one degree per second 
T 

with e^ (0) = (1,0, 0,0). The values of 3c(l) are computed for 

(1 - 1 ) 

p = 2, 3,... 10 and T=2 V ' , i = 1,2.. *10, using double precision 
(14 hexadecimal digits) . The theoretical values of eXl) are determined 
from equation (11-13) 
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e^U) 


Cos (V3 tt/360) 

e 2 a) 


(l/yj) Sin(./3 tt/360) 

e 3 a) 


(1/./3) Sin(i/3 tt/360) 

e 4 (l) 


(1/V3) Sin 0/3 ir/360) 


The computed values of e_(l) and the theoretical values of e/1) are com- 
pared so as to obtain the actual computation error* Some of the re- 
sulting actual computational error norms are plotted in Fig* 4 as a 
function of the time increment T for p = 2,3,4 and 7. Notice that be- 
tween T = 1 and T = .001, the computational error is dominated by the 
truncation error for p = 2 and is dominated by the roundoff error for 
p = 7. It is observed that minimum computational error norm occurs at 
T = 2~^ second, T = , T = 2”^, and T = 2 ^ for p = 3, p = 4, p=5, 

and p = 6, respectively. 

Theoretical Computational Error Norm 

The norm of the truncation error is obtained from equation (III-19) 
and the norm of the roundoff error is obtained from equation (III-131) • 
The norm of the theoretical computional error is computed by adding the 
truncation and roundoff error norm. Some of the results are depicted 
in Figures 5 through 7. 

Fig. 5 shows the theoretical truncation error norm and theoretical 
roundoff error norm as a function of the time increment for p - 2, 4 
and 6. Again, it illustrates all the characteristics described in 
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Figure 4. Actual computational error norms vs. time increment T. 



Error Norm 
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Example 1. It may be noted that the roundoff error is not as sensitive 

to the order of the finite series, p, as the truncation error is. 

Fig. 6 shows the theoretical computational error norm as a function 

of the time increment T. It is observed that the optimal T for p » 3, 

—6 *-4 

p = 4, p = 5, and p - 6 is 2~ second, 2 ~ second, second, and 1 
second , respectively, 

! 

In Fig, 7, theoretical curves for the computational error norm are 
compared with the actual curve for p - 2,3,4 and 7. It may be seen that 
using higher order finite series with a larger time increment will re- 
duce the speed for the computer in calculating as well as decrease 
the computational error. 



V. CONCLUSIONS AND RECOMMENDATIONS 


A numerical integration scheme (the finite series method) for 
solving the four parameter vector differential equation is derived 
and investigated in this report. The results obtained can be applied 
to a large class of numerical integration schemes, since this class 
can be shown to be equivalent to the finite series approximation 
method . 

Studies show that there are two types of computational errors in 
computing the numerical solutions to the four parameter vector 
differential equation using a digital computer. These are truncation 
error and roundoff error. Truncation error is caused by the 
approximate nature of the numerical integration scheme. Roundoff 
error is due to the fact that all numbers are represented by a finite 
number of digits in a computer. ' 

Bounds for the truncation errors and roundoff errors generated 
by the computer in computing the four parameters using the finite 
series method are derived. The results show that the truncation 
error norm can be expressed as a function of the initial conditions 
of the four parameters, the magnitude of the angular rotations and 
the number of terms used in the numerical approximation of the state 
transition matrix. The results also illustrate that the roundoff 
error norm can be expressed as a function of the initial conditions 
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of the four parameters, the magnitude of the angular rotation and the 
local roundoff error. The local roundoff error in turn can be 
expressed as a function of the number system and number of digits 
employed by the digital computer and the number of terms used -in the 
numerical approximation td the state transmission matrix. Study 
results show that the error norm developed is useful in the determina- 
tion of an optimal integration step size for the four parameter algorithm, 
and the computer sizing requirement for a particular mission. 

It should be emphasized that the computational error norm derived 
in this analysis is an upper bound on the error generated by the 
digital computer in computing the four parameters using the finite 
series method. The actual errors that would be observed might therefore 
be and are shown to be considerably less than the error analytically 
determined by this method. Nevertheless, this technique does provide 
a means of obtaining the limit that can be placed on the errors caused 
by the computational process using a digital computer. In order to 
obtain a more realistic bound on the roundoff error, a statistical 
approach should be investigated. 
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APPENDIX A 


SOLUTION OP THE POUR PARAMETERS USING THE 
PEANO -BAKER METHOD OF 
SUCCESSIVE APPROXIMATION 19] 


Consider e = A(t) e(t) 


(A-l) 


Integrating (A-l) gives : 


e.(t) = e.(t ) + f t A(t) £(t)dT 


(A-2) 


Equation (A-2) may be solved by an iterative scheme called the Peano- 
Eaker method of successive approximations which involves repeated sub- 
stitution of e_(t) from the left member of (A-2) into the integral. 

£ 

1 st iteration: _e(t) = e/t ) + / A(t) e.(t ) dr 


= {I + f A(T)dT}e.(t 0 ) 
C o 


2 n xteration: e(t) = e(t ) + / A(t) {I + / A(r)dT}e(t )dr 
— — o J t t ~o 

u o o 

t t t 

- {1 + f A(x)dr + / A(x) [/ A(T)dT]dT}e^(t Q ) 
to t D t Q 


Thus an infinite series can be obtained. If the elements of the 
matrix A(t) remain bounded in the range from 0- to t, it may be shown that 
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the infinite series converges to the solution e_(t) . 

If the elements of A(t) are of the form let A(t) = Bt 

where 


t. 



-a 


3 


-a 


2 


-a. 


B 



a 3 

a 2 


0 

a l 


~ a l 


0 


a 2 


“ a 3 


a l 


a 


2 


a 


3 


0 


then 

t t t 

e(t) = {I + / A(x)dx + / A(x)l/. A(x)dx]dT + — } e(t ) 

to fc o “ ° 


“ {I + / t (Bx)dt + J t (Bx)dx]dx + ...} e.(t Q ) 


t t t 

=* {I + B / Tdx + B^ / t f / TdxJdT + e(t ) 

v t t't — o 

o o o 


Let t - 0 
o 


1 2 v?*** B n t^ n 

eCt) - {I + 7 Bt 2 + + . . . ■ ■ - 

2 2 x 4 2 n (nl) 


+ . e(t Q ) 


; 1 _.2 . B 1 2 t 4 


B 3 t« 


bV 


bV° 


{I + "2 Bt + — + -^ + ? _ + 5 ^ T+ ... } e(t> 


Since B is a skew-symmetric matrix, the following identities can be 


obtained 
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B 


,1.2, 2 2 2. _ 

v»p a 2 


- c 2 I 


where c 2 = (a 2 + a^ + a 2 ) 


c|) 2 = k 2 • c^) 2 


B = 


tC 2 B 


4 4 

B = cl 


In general 


n 


n-1 


c = (-D 2 c n 1 B 


for n odd 


n 


= (-1) 2 c n I 


for n even 


Thus 


2 c 2 .(|i) 3 c 4 (f) 5 

• (t) -~if — + — if — 




+ I[1 - 


2 t 2 2 4 t 2 4 

c 2 e~) c cf-r 


2 ! 


4! 


- ...]} e(t ) 


3,1?, 3 5,t!,5 


. B . 2 c J (±-) 


...] 


2 

+ l{Cos (^~ )]} ^(t Q ) 

2 2 

=' {f Sin(~~) + I [Cos (—■“—) ] } e(0 
2 2 

e(t) = {| S±nC~ ) + I[Cos(^|- )]>e(t D ) 


(A-3) 


Equation (A-3) is the exact solution for the elements of the e vector 

clt 

if the elements of - A in equation (A-l) are of the form (“) over the time 


interval 0 to t. 



APPENDIX B 


VECTOR AND MATRIX NORMS 

The norm of an N-vector x is a real, non-negative number, 
denoted by | jxj |, which gives an assessment of the size of the vector. 
This norm satisfies the following properties 


| |x| | > 0 if x # £ (B-l) 

j |x| | = 0 if x= 0 (B-2) 

I i k *l I = ! k l I lil I where k is a scalar (B-3) 

! 1^1 I 1 11*1 I + I !xl I ( B "4) 

Prom inequality (B-4) , the following inequality is deduced 

I Ix-X| I > | |i| | - I |X| | ( B "5) 

i I-Xl I > I Ixl I - I lil I < B -6) 


The most commonly used vector norms are defined by 
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(2) I|x| 1 2 " E I } and 

i=l 

^ I 1*1 L = max lail 

i 

Similarly, the norm of an (NxN) -matrix A is a real, non-negative 
number, denoted by | J a| | , which satisfies 


| | A J | > 0 if A 4 [0] - (B-7) 

! I A | | = 0 if A = [0] (B- 8 ) 

! 1^1 I = l k ! I |A| | where k is a scalar - (B-9) 

I l A+B l I 1 \ l A l I + I l B [ [ (B— 10) 

ll*sll <IWI 11*11 CB- ii) 

IMI I IWI l|B|| (B-12) 


The matrix norms corresponding to the 1,2 and <»— vector norms are, 
respectively: 

(D [Ml-j^max l |a i3 j (B-13) 

j i=l 


(2) | | A J j 2 ~ (maximum eigenvalue of A H A)l/2 


. (B-14) 



where 
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denotes the complex conjugate transpose of A } 


and 


(3) | |A| | ro = max £ \a ± A 
i j=l 3 


(B-15) 



APPENDIX C 


To prove that for proportional angular rates 


<!>(m)E(k) = E(k)$(m), for all positive integers k & m (C-l) 


first consider ®(m)$(k). 

Prom equation (11-10) 

,, ^ v [a(m)K]^ 

®(m) = I TT 

i=0 1 * 

S(k) = l UC^K ] 1 
i=0 i! 


, it can be shown that 


and 


(C-2) 


(C-3) 


where K is a constant matrix defined by equation (II-4) and a(k) and 
a(m) are scalar functions. Thus 


$(m)$(k) ={ l 

i=0 


[a(m)K] 



[a(k)K] 1 

i! 


(C-4) 


Since a(m) and a(k) are scalar functions and K is a constant matrix, then 


{ l [a(nt)Kl 1 1 . [a(k)K1^ = [a(k)K] 1 { £ [a(m)K]S ( C -5) 

• i=0 i I i ! i ! i=0 i ! 
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Therefore 


°° • p , 

$(m)$(k) = y { Ia(m)K] x j { V [aCkjK] 1 } 
i=0 ii i=0 il 

= f ' { [a(k)K] i \ r y [a(m)K]^ -i 

• /s •{ » i 4 t J 

r=0 x * i=0 x * 

= $(k)*(m) (C-6) 

Next consider $(m)$(k). For proportional angular rate, 

$(m)$(k) is defined as 


^ K ' ^ K 


^(m)$(k) = e z 

(C-7) 

Since = [^4] [^k] 

(C-8) 

a (m) a(k) a(k) a(m) 

- . * . O ^ 0 0 K rj JX 



then <Km)$(k) = e 2 e 2 =e 2 e 2 = $(k)$(m) (C-9) 


Now consider $(m)E(k). From equation (III-25) , $(m)E(k) can be 
represented as 


$(m)E(k) = *(m)[$(k). - $(k)J = $(m)$(k) - $(m)$(k) 


(C-10) 



88 


Substituting equations (C-6) and (C-9) into equation (C-10) yields 


$(m)ECk) = [0(k) - S(k>] $(m) = E(k)$(m) 


(C-ll) 



APPENDIX D 


DIGITAL COMPUTER PROGRAM FOR SOLUTION 
OF EXAMPLE 1 IN CHAPTER IV 


DIMENSION TX(2) , A(2,2) ,B(2,2) ,C(2,2)D(2,2)PHI(2,2) 
1,F(20,T,DT,X(2),X0(2) 

DOUBLE PRECISION R 

TX (1) =2 . *DEXP (-1 . DO ) -DEXP (-2 . DO)+DEXP (-1 . DO) -DEXP'(-2 . DO) 
TX(2)=2 . * (DEXP (-2 . DO) -DEXP (-1 .DO) )+2 . *DEXP (-2 . DO) -DEXP (-1 .DO) 
WRITE (6,22) TX(1) ,TX(2) 

22 FORMAT (IX , 6HTX (l')= , E16 . 8 , 5X , 6HTX (2) = , E16 . 8) 

DO 100 K=5 , 11 

F(2)=2.EO 

F(3)=6.EO 

F(4)=2.4E1 

F(5)=1.2E2 

F(6)=7.2E2 

F(7)=5.04E3 

F(8)=4.032E4 

F(9)=3.6288E5 

F(10)=3.6288E6 

F(11)=3.99168E7 

DO 100 M=l,19 

XO(l)=l, 

XO(2)=l. 

T=.0005 
DT=2 . ** (1-M) 

TST0P=1 . 

A(1,1)=0. 

A(1,2)=1.*DT 

A(2,1)=2.*DT 

A(2,2)=3.*DT 

DO 1 1=1,2 

DO 1 J-1,2 

D(I,J)=0.0 

D(I,I)=1.0 

PHI(I,J)=D(I,J)+A(I,J) 

1 C(I, J)=A(I, J) 

N=2 

DO 3 L=2,K 
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DO 2 1=1,2 
DO 2 J=l,2 

2 B(I,J)=C(I,J) 

CALL MATMUL (A,B,N,C) 

DO 3 1=1,2 
DO 3 J=l,2 

3 PHI (I , J) =PRI (I , J )+C (I , J) /F (J) 

5 T=T+DT 

DO 7 1=1,2 
X(I)=0. 

DO 7 J-1,2 

7 X (I) =PHI (I , J) *XO (J)+X (I) 

DO 8 1=1,2 

8 XO(I)=X(I) 

IF(T.LT.TSTOP) GO TO 5 

R= ABS(X(1)-TX(L))+ ABS (X(2)-TX(2)) 

WRITE (6, 23) DT,K 

23 FORMAT (IX , 3HDT= , F6 . 3 , 2HK= , 12 ) 

DO 6 1=1,2 

6 WRITE (6,21) T,I,X(I) 

21 FORMAT (2X, 2HT= , F6 . 3 , 7X, 2HX ( , 12 , 2H) = , E16 . 8) 

WRITE (6 ,42) R 

42 FORMAT (10X,2HR=,D23.16 ,//) 

100 CONTINUE 
STOP 
END 


SUBROUTINE MALMUL (A,B,N,C) 

DIMENSION A(2,2),B(2,2),C(2,2) 

DO 1 1=1, N 

DO 1 J=1,N 

C(I,J)=0.0 

DO 1 K=1,N 

1 C(I,J)=C(I,J)+A(I,K)*B(K,J) 

RETURN 

END 



APPENDIX E 


DIGITAL COMPUTER PROGRAM TO COMPUTE 
THE FOUR PARAMETERS USING THE 
FINITE SERIES METHOD 


DOUBLE PRECISION A(4,4) ,B (4,4) ,C(4,4) ,OMEGA(4,4) , 
IE (4 ) , EO (4 ) , PHI ( 4 , 4 ) PHIXD , PHIYD ,PHIZD , T , D (4 , 4) , DT , 
10MEGO(4,4),E1,E2,AA,CC,EE,R,F(20) 

F(2)=2.EO 

F(3)=6.E0 

F(4)=2.4E1 

F(5)=1.2E2 

F(6)=7.2E2 

F(7)=5.04E3 

F(8)=4.032E4 

F(9)=3.6288E5 

F(10)=3.6288E6 

F (11)=3 . 99168E7 

AA=3 . * (3.14159265/180.) **2. 

CC=DSQRT(AA) 

EE=CC/2. 

E1=DC0S (EE) 

E2=DSIN (EE) /DSQRT (3 . DO) 

WRITE (6.41) E1,E2 

41 FORMAT (10X, 3HE1= ,D23 . 16 , 2X, 3HE2= ,D23 . 16) 

PHIXD=3 . 14159265/180 . 

PHIYD= 3.14159262/180. 

PHIZD=3 . 14159265/180 . 

OMEGA (1,1)=0.0 
0MEGA(1,2)=-PHIZD 
OMEGA (1 , 3) =-PHIYD 
OMEGA (1,4) =-PHIXD 
OMEGA (2,1) =PHIZD 
OMEGA (2, 2) =0.0 
OMEGA (2 , 3) =-PHIXD 
OMEGA(2,4)= PHIYD 
0MEGA(3,1)= PHIYD 
OMEGA (32,)= PHIXD 
0MEGA(3,3)=0.0 


91 



noooo ono ooo 


92 


0MEGA(3,4)=-PHIZD 
OMEGA (4,1) =PHIXD 
OMEGA (4 , 2 )=— PHIYD 
OMEGA (4,3)= PHIZD 
OMEGA (4 ,4) =0.0 
OMEGO (1 , 1)=0 . 0 
OMEGO (1 , 2) =0 . 0 
OMEGO (1,3)=0.0 
OMEGO (1,4)=0.0 
OMEGO (2,1) =0.0 
OMEGO (2, 2) =0.0 
OMEGO (2, 3) =0.0 
OMEGO(2,4)=0.0 
OMEGO(3,1)=0.0 
OMEGO (3,2)=0.0 
OMEGO(3,3)=O.0 
OMEGO(3,4)=O.0 
OMEGO (4 , 1)=0 . 0 
OMEGO (4,2)=0.0 
OMEGO (4 , 3)=0 . 0 
OMEGO(4,4)=0.0 
DO 100 K=2,10 
DO 100 M=1 , 10 
T= . 0005 
DT=2.**(1-M) 

TSTOP=l . 

FO (1)=1.0 
E0(2)=0. 

E0(3)=0. 

E0(4)=0. 

DO 1 1=1,4 
DO 1 J-1,4 

A (I , J) =DT* (OMEGA (I,J) )/2. 

CALCULATION OF STATE TRANSITION MATRIX PHI 

D(I, J)=0 .0 
D(I,I)=1.0 


CALCULATION OF THE FIRST TWO TERMS OF THE STATE 
TRANSITION MATRIX 
PHI(I, J)=D (I, J)+A(I, J) 

1 C(I, J)=A(I, J) 

CALC. OF PHI FOR P>2 
N IS THE ORDER OF THE SYSTEM 

(K=1)=NUMBER OF TERMS USED IN THE INFINITE SERIES 
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N=4 

DO 3 L=2,K 
DO 2 1-1,4 
DO 2 J=l,4 

2 B(I,J)=C(I,J) 

CALL MATMUL (A,B,N,C) 

DO 3 1=1,4 
DO 3 J=l,4 

3 PHI(I,J)=PHI(I,J)+C(I,J)/F(L) 

5 T=T+DT 

DO 7 1=1,4 
E(I)=0. 

DO 7 J=l,4 

7 E (I) =PHI (I , J ) *EQ ( J )+E (I ) 

DO 8 1=1, N 

8 EO(I)=E(l) 

IF(T.LT.TSTOP) GO TO 5 

R=DABS (E(1)-E1)+DABS (E2)+DABS (E(3)-E2) 
1+DABS (E(4)-E2) 

WRITE (6, 22) DT,K 

22 FORMAT (IX , 3HDT= , F6 . 3 , 2HK= ,12) 

DO 6 1=1,4 

6 WRITE (6,21) T,I ,E(I) 

21 FORMAT (2X,2HT=,F6 .3, 7X,2HE(,I2 ,2H)=,D23.16) 
WRITE (6, 42) R 

42 FORMAT (10X,2HR=,D23.16,//) 

100 CONTINUE 
STOP 
END 


SUBROUTINE MATMUL (A,B,N,C) 

DOUBLE PRECISION A(4,4) ,B(4,4) ,C(4,4) 
C CALCULATE C(I,J) COEFFICIENTS 
DO 1 1=1, N 
DO 1 J=1,N 
C(I,J)=0.0 
DO 1 K=1,N 

1C (I, J)=C(I, J)+A(I,K) *B(K, J) 

RETURN 

END 



